This is the analysis of all the salinity and nutrient experiments conducted on the lanai in St. John 616 from September 2021 to October 2022. These experiments incorporated four paired salinity and nutrient treatments with three temperatures. Each of the first four runs produced an n = 2 and was repeated initially 8 times for a total of n = 16. Data gaps were identified and filled in February, April, and October 2022. This output reflects all data totaling five treatments for Ulva and six treatments for Hypnea.

Packages loaded:

library(lme4)
library(lmerTest)
library(afex)
library(effects)
library(car)
library(MuMIn)
library (dplyr)
library(emmeans)
library(DHARMa)
library(performance)
library(patchwork)
library(rstatix)
#for plots and tables
library(ggplot2)
library(ggpubr)
library(forcats)
library(RColorBrewer)
library(tidyverse)
library(sjPlot)
library(sjmisc)
library(mmtable2)
library(gt)
library(purrr)
library(stringr)
library(tidyr)

Load this file for normalized to quantum efficiency of photosynthesis per same as above

all_runs_photosyn_data <- read.csv("../data_input/hyp_ulva_all_runs_ek_alpha_normalized.csv")

Assign numerous variables as factors

all_runs_photosyn_data$Run <- as.factor(all_runs_photosyn_data$Run)
all_runs_photosyn_data$Temperature <- as.factor(all_runs_photosyn_data$Temp...C.)
all_runs_photosyn_data$Treatment <- as.factor(as.character(all_runs_photosyn_data$Treatment))
all_runs_photosyn_data$deltaNPQ <- as.factor(all_runs_photosyn_data$deltaNPQ)

Subset data to get only Ulva data for day 9 Ek– also remove the treatment 2.5 that was problematic

ulva <- subset(all_runs_photosyn_data, Species == "ul" & RLC.Day == 9 & Treatment != 2.5)
ulva$treatment_graph[ulva$Treatment == 0] <- "1) 35ppt/0.5umol"
ulva$treatment_graph[ulva$Treatment == 1] <- "2) 35ppt/14umol" 
ulva$treatment_graph[ulva$Treatment == 2] <- "3) 28ppt/27umol" 
ulva$treatment_graph[ulva$Treatment == 3] <- "5) 18ppt/53umol" 
ulva$treatment_graph[ulva$Treatment == 4] <- "6) 11ppt/80umol"
ulva$treatment_graph[ulva$Treatment == 2.5] <- "4) 28ppt/53umol"

Ulva Ek________________________________________________________________________________________

hist(ulva$ek.1, main = paste("Ulva lactuca"), col = "olivedrab3", labels = TRUE)

ulva %>% ggplot(aes(ek.1)) +
  geom_histogram(binwidth=5, fill = "#5BB300", color = "black", size = 0.25, alpha = 0.85) +
  theme_bw()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.

Ulva Ek model

ulva_ek_model <- lmer(formula = ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva)

Make residual plots of the data for ulva ek

hist(resid(ulva_ek_model))

plot(resid(ulva_ek_model) ~ fitted(ulva_ek_model))

qqnorm(resid(ulva_ek_model))
qqline(resid(ulva_ek_model))

Check the performance of the model, get r2, make a table, plot effects, summary

performance::check_model(ulva_ek_model)
## Variable `Component` is not in your data frame :/

r.squaredGLMM(ulva_ek_model)
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
##            R2m       R2c
## [1,] 0.6248358 0.7151076
summary(ulva_ek_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +  
##     (1 | RLC.Order)
##    Data: ulva
## 
## REML criterion at convergence: 2145.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.75984 -0.58531  0.07584  0.56071  2.65765 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Plant.ID  (Intercept)  76.53    8.748  
##  Run       (Intercept)  43.24    6.576  
##  RLC.Order (Intercept)  18.86    4.343  
##  Residual              437.52   20.917  
## Number of obs: 240, groups:  Plant.ID, 95; Run, 7; RLC.Order, 6
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)     34.248      6.522   8.394   5.251 0.000661 ***
## Treatment1      39.125      7.213   6.771   5.425 0.001097 ** 
## Treatment2      46.727      7.213   6.771   6.479 0.000392 ***
## Treatment3      79.171      7.213   6.771  10.977 1.47e-05 ***
## Treatment4      84.260      7.213   6.771  11.682 9.80e-06 ***
## Temperature27  -10.542      4.842  32.140  -2.177 0.036904 *  
## Temperature30  -10.619      4.435  63.820  -2.394 0.019600 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Trtmn2 Trtmn3 Trtmn4 Tmpr27
## Treatment1  -0.688                                   
## Treatment2  -0.688  0.825                            
## Treatment3  -0.688  0.825  0.825                     
## Treatment4  -0.688  0.825  0.825  0.825              
## Temperatr27 -0.355  0.002  0.002  0.002  0.002       
## Temperatr30 -0.344  0.000  0.000  0.000  0.000  0.475
tab_model(ulva_ek_model)
  ek.1
Predictors Estimates CI p
(Intercept) 34.25 21.40 – 47.10 <0.001
Treatment [1] 39.12 24.91 – 53.34 <0.001
Treatment [2] 46.73 32.52 – 60.94 <0.001
Treatment [3] 79.17 64.96 – 93.38 <0.001
Treatment [4] 84.26 70.05 – 98.47 <0.001
Temperature [27] -10.54 -20.08 – -1.00 0.030
Temperature [30] -10.62 -19.36 – -1.88 0.017
Random Effects
σ2 437.52
τ00 Plant.ID 76.53
τ00 Run 43.24
τ00 RLC.Order 18.86
ICC 0.24
N Run 7
N Plant.ID 95
N RLC.Order 6
Observations 240
Marginal R2 / Conditional R2 0.625 / 0.715 
plot(allEffects(ulva_ek_model))

Construct null model to perform likelihood ratio test REML must be FALSE

ulva_ek_treatment_null <- lmer(formula = ek.1 ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
## boundary (singular) fit: see help('isSingular')
ulva_ek_model2 <- lmer(formula = ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
anova(ulva_ek_treatment_null, ulva_ek_model2)
## Data: ulva
## Models:
## ulva_ek_treatment_null: ek.1 ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## ulva_ek_model2: ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                        npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## ulva_ek_treatment_null    7 2332.8 2357.1 -1159.4   2318.8                     
## ulva_ek_model2           11 2200.1 2238.3 -1089.0   2178.1 140.69  4  < 2.2e-16
##                           
## ulva_ek_treatment_null    
## ulva_ek_model2         ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ulva_ek_temperature_null <- lmer(formula = ek.1 ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
ulva_ek_model3 <- lmer(formula = ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
anova(ulva_ek_temperature_null, ulva_ek_model3)
## Data: ulva
## Models:
## ulva_ek_temperature_null: ek.1 ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## ulva_ek_model3: ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                          npar    AIC    BIC  logLik deviance  Chisq Df
## ulva_ek_temperature_null    9 2202.8 2234.2 -1092.4   2184.8          
## ulva_ek_model3             11 2200.1 2238.3 -1089.0   2178.1 6.7989  2
##                          Pr(>Chisq)  
## ulva_ek_temperature_null             
## ulva_ek_model3              0.03339 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plot raw data

ulva %>% ggplot(aes(treatment_graph, ek.1)) + 
  geom_boxplot(size=0.5) + 
  geom_point(alpha = 0.75, size = 3, aes(color = Temperature), show.legend = TRUE) + 
  labs(x="salinity/nitrate", y= "Day 9 Ek (μmols photons m-2 s-1)", title= "A", subtitle = "Ulva lactuca") + 
  scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "18ppt/53umolN", "11ppt/80umolN")) + 
  ylim(-1, 250) + stat_mean() +
  scale_color_manual(values = c("#295102", "#7CB950", "#BDE269")) +
  geom_hline(yintercept=0, color = "red", size = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
        plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))

Make linear regression with growth add growth rate from other dataset to this one and subset by species

growth_rate <- read.csv("/Users/Angela/src/work/limu/algal_growth_photosynthesis/data_input/all_runs_growth_011723.csv")
growth_rate$Species <- as.factor(growth_rate$Species)
growth_rate$treatment <- as.factor(growth_rate$treatment)
growth_rate$lunar.phase <- as.factor(growth_rate$lunar.phase)

Make a new column for weight change (difference final from initial) add growth column to dataset

growth_rate$growth_rate_percent <- (growth_rate$final.weight - growth_rate$Initial.weight) / growth_rate$Initial.weight * 100
gr_ulva <- subset(growth_rate, Species == "Ul" & treatment != 2.5)

ulva$growth_rate <- round((gr_ulva$final.weight - gr_ulva$Initial.weight) / gr_ulva$Initial.weight * 100, digits = 2)
ulva$lunar.phase <- (gr_ulva$lunar.phase)

Plot linear regression between ek and growth

ulva_growth_ek_graph <- ggplot(ulva, aes(x=ek.1, y=growth_rate)) + 
  geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = Treatment)) + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Ulva lactuca Ek vs Growth Rate", x = "Ek (μmols photons m-2 s-1)", 
       y = "growth rate (%)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor()
ulva_growth_ek_graph
## `geom_smooth()` using formula = 'y ~ x'

Summarize the means for Ek

ulva %>% group_by(Treatment) %>% summarise_at(vars(ek.1), list(mean = mean))
## # A tibble: 5 × 2
##   Treatment  mean
##   <fct>     <dbl>
## 1 0          27.2
## 2 1          65.2
## 3 2          72.8
## 4 3         105. 
## 5 4         110.

Hypnea Ek____________________________________________________________________ Use Day 9 for final analysis hm1-1 on 10/12/22 and hm1-2 on 4/29/22 causing issues of influential observations hm1-2 for both rETRmax and Ek – leaving them in dataset because no good reason to believe not good data There was no D9 RLC for hm6-4 on 11/12/21 but had to remove hm6-4 from 10/9/21 below to match growth data

hypnea <- subset(all_runs_photosyn_data, Species == "hm" & RLC.Day == 9 & uid != "2021-10-09_hm6-4")
hypnea$treatment_graph[hypnea$Treatment == 0] <- "1) 35ppt/0.5umol"
hypnea$treatment_graph[hypnea$Treatment == 1] <- "2) 35ppt/14umol" 
hypnea$treatment_graph[hypnea$Treatment == 2] <- "3) 28ppt/27umol" 
hypnea$treatment_graph[hypnea$Treatment == 3] <- "5) 18ppt/53umol" 
hypnea$treatment_graph[hypnea$Treatment == 4] <- "6) 11ppt/80umol"
hypnea$treatment_graph[hypnea$Treatment == 2.5] <- "4) 28ppt/53umol"

Make a histogram for hypnea

hist(hypnea$ek.1, main = paste("Hypnea musciformis -- Ek"), col = "maroon2", labels = TRUE)

hypnea %>% ggplot(aes(ek.1)) +
  geom_histogram(binwidth=5, fill = "#7D0033", color = "black", size = 0.25, alpha = 0.85) +
  theme_bw()

Run model for Hypnea and Ek

hyp_ek_model <- lmer(formula = ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID), data = hypnea)

Plot residuals

hist(resid(ulva_ek_model))

plot(resid(hyp_ek_model) ~ fitted(hyp_ek_model))

qqnorm(resid(hyp_ek_model))
qqline(resid(hyp_ek_model))

Check the performance of the model, get r2, make a table and plot for model

performance::check_model(hyp_ek_model)
## Variable `Component` is not in your data frame :/

r.squaredGLMM(hyp_ek_model)
##            R2m       R2c
## [1,] 0.3763202 0.5670999
summary(hyp_ek_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID)
##    Data: hypnea
## 
## REML criterion at convergence: 2625.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0609 -0.6195 -0.0476  0.4318  3.7844 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Plant.ID (Intercept) 112.8    10.62   
##  Run      (Intercept) 127.8    11.31   
##  Residual             546.1    23.37   
## Number of obs: 286, groups:  Plant.ID, 143; Run, 8
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)     76.847      9.143   6.026   8.405 0.000151 ***
## Treatment1      18.182     10.840   6.057   1.677 0.144009    
## Treatment2      17.509     10.840   6.057   1.615 0.156904    
## Treatment2.5    42.029     14.806   4.602   2.839 0.039907 *  
## Treatment3      17.303     10.840   6.057   1.596 0.161073    
## Treatment4      24.523     10.864   6.113   2.257 0.063977 .  
## Temperature27  -35.829      4.237 102.763  -8.456 1.98e-13 ***
## Temperature30  -39.418      4.224 106.633  -9.331 1.71e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Trtmn2 Trt2.5 Trtmn3 Trtmn4 Tmpr27
## Treatment1  -0.784                                          
## Treatment2  -0.784  0.903                                   
## Treatmnt2.5 -0.574  0.484  0.484                            
## Treatment3  -0.784  0.903  0.903  0.484                     
## Treatment4  -0.783  0.901  0.901  0.483  0.901              
## Temperatr27 -0.231  0.002  0.002  0.000  0.002  0.002       
## Temperatr30 -0.231  0.000  0.000  0.000  0.000  0.005  0.496
tab_model(hyp_ek_model)
  ek.1
Predictors Estimates CI p
(Intercept) 76.85 58.85 – 94.85 <0.001
Treatment [1] 18.18 -3.16 – 39.52 0.095
Treatment [2] 17.51 -3.83 – 38.85 0.107
Treatment [2.5] 42.03 12.88 – 71.18 0.005
Treatment [3] 17.30 -4.04 – 38.64 0.112
Treatment [4] 24.52 3.14 – 45.91 0.025
Temperature [27] -35.83 -44.17 – -27.49 <0.001
Temperature [30] -39.42 -47.73 – -31.10 <0.001
Random Effects
σ2 546.12
τ00 Plant.ID 112.83
τ00 Run 127.85
ICC 0.31
N Run 8
N Plant.ID 143
Observations 286
Marginal R2 / Conditional R2 0.376 / 0.567 
plot(allEffects(hyp_ek_model))

Construct null model to perform likelihood ratio test REML must be FALSE

hyp_ek_treatment_null <- lmer(formula = ek.1 ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
hyp_ek_model2 <- lmer(formula = ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
anova(hyp_ek_treatment_null, hyp_ek_model2)
## Data: hypnea
## Models:
## hyp_ek_treatment_null: ek.1 ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## hyp_ek_model2: ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                       npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## hyp_ek_treatment_null    7 2690.8 2716.3 -1338.4   2676.8                     
## hyp_ek_model2           12 2690.0 2733.8 -1333.0   2666.0 10.789  5    0.05574
##                        
## hyp_ek_treatment_null  
## hyp_ek_model2         .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hyp_ek_temp_null <- lmer(formula = ek.1 ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
hyp_ek_model3 <- lmer(formula = ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
anova(hyp_ek_temp_null, hyp_ek_model3)
## Data: hypnea
## Models:
## hyp_ek_temp_null: ek.1 ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## hyp_ek_model3: ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                  npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## hyp_ek_temp_null   10 2740.4 2777.0 -1360.2   2720.4                         
## hyp_ek_model3      12 2690.0 2733.8 -1333.0   2666.0 54.454  2  1.498e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Use ggplot for the data in a boxplot

hypnea %>% ggplot(aes(treatment_graph, ek.1)) + 
  geom_boxplot(size=0.5) + 
  geom_point(alpha = 0.75, size = 3, aes(color = Temperature), show.legend = TRUE) + 
  labs(x="salinity/nitrate", y= "Ek (μmols photons m-2 s-1)", title= "B", subtitle = "Hypnea musciformis") + 
  scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) + 
  ylim(-1, 250) + stat_mean() + 
  scale_color_manual(values = c("#9C0627", "#BB589F", "#F4B4E2")) +
  geom_hline(yintercept=0, color = "red", size = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
        plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))

For Hypnea remove hm6-4 on 11/12 that had no d9 RLC (final weight 0.1017) and hm6-4 on 10/9/21 because it was white and also looked dead

gr_hypnea <- subset(growth_rate, Species == "Hm" & final.weight != 0.1017 & growth_rate_percent > -87.96837)
hypnea$growth_rate <- round((gr_hypnea$final.weight - gr_hypnea$Initial.weight) / gr_hypnea$Initial.weight * 100, digits = 2)

Plot a regression between the photosynthetic independent variables of interest and growth rate

hypnea_growth_ek_graph <- ggplot(hypnea, aes(x=ek.1, y=growth_rate)) + 
  geom_point(alpha = 0.5, size = 3, aes(color = Treatment), show.legend = TRUE) + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Hypnea musciformis Ek vs Growth Rate", x = "Ek (μmols photons m-2 s-1)", 
       y = "growth rate (%)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor(label.y = 150)
hypnea_growth_ek_graph
## `geom_smooth()` using formula = 'y ~ x'

Summarize the means

hypnea %>% group_by(Treatment) %>% summarise_at(vars(ek.1), list(mean = mean))
## # A tibble: 6 × 2
##   Treatment  mean
##   <fct>     <dbl>
## 1 0          51.8
## 2 1          71.4
## 3 2          70.7
## 4 2.5        93.8
## 5 3          70.5
## 6 4          78.3
hypnea %>% group_by(Run) %>% summarise_at(vars(ek.1), list(mean = mean))
## # A tibble: 8 × 2
##   Run    mean
##   <fct> <dbl>
## 1 1      90.1
## 2 2      79.1
## 3 3      60.8
## 4 3.5    64.5
## 5 4      58.7
## 6 7      93.8
## 7 8      52.2
## 8 9      51.3